
#include <complex>
#include <fstream>
#include <iostream>
#include <string.h>
#include <iomanip>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include<time.h>
#include<iomanip>
#include <sstream>


using namespace std;


#define pi 3.1415926535898
#define CONFINED

template<typename T>inline int anint(T x)
	{  return(int(floor(x+0.5)));  }

struct PART {
   double x;
   double y;
   double z;
   double sigma;
   double energia;
};

struct BOX{
    double lx;
    double ly;
    double lz;
    int npart;
};


int main(void)
{
    void SetRandomPart(BOX bx, PART *part, double phi, double poly, double ratio);
    int npart=1000;
    PART *part = new PART[npart];
    double poly=0.;
    double ratio=1.0;
    double phi=0.2;
    BOX box;
    box.npart=npart;
    srand (time(NULL));
    SetRandomPart(box, part, phi, poly, ratio);

    ofstream MiArchivo ("posicion.dat",ios::out);
	if (MiArchivo.is_open()) {
	    for (int i=0;i<box.npart;i++)
	    {
                MiArchivo<<part[i].x<<" "<<part[i].y<<" "<<part[i].z<<" "<<part[i].sigma<<endl;
	    }
	}
	MiArchivo.close();

    delete []part;
    return 0;
}


void SetRandomPart(BOX bx, PART *part, double phi, double poly, double ratio)
{
    double distancia2(BOX bx, PART *part, int i, int j);
    int nok=0;
    bool good;
    double r2;
    double dmin;

    //double poly=0.36;
    double Vpart=0.0;
    double vol;
    double u1,u2;
    for(int i=0;i< bx.npart;i++)
    {
        u1=double(rand())/RAND_MAX;
        u2=double(rand())/RAND_MAX;
        part[i].sigma=exp(poly*sqrt(-2.*log(u1))*cos(2.*pi*u2));
        Vpart+=1./6.*pi*pow(part[i].sigma,3.0);
    }
    vol=Vpart/phi;

    bx.lx=pow(vol/ratio,1./3.);
    bx.ly=bx.lx;
    bx.lz=ratio*bx.lx;
    cout<<bx.lx<<"  "<<bx.ly<<"  "<<bx.lz<<endl;
    dmin=bx.lx;

	for(int i=0;i< bx.npart;i++)
	{
        part[i].x=bx.lx*(double(rand())/RAND_MAX);
        part[i].y=bx.ly*(double(rand())/RAND_MAX);
        //part[i].z=bx.lz*(double(rand())/RAND_MAX);
        part[i].z=(bx.lz-2.*part[i].sigma)*(double(rand())/RAND_MAX)+part[i].sigma;
        //part[i].z=(bx.lz-2.*part[i].sigma)*(double(rand())/RAND_MAX)+part[i].sigma;

		///Comprobamos solapamiento con otras particulas
		good=true;
		double daux=1.e7;
		for(int j=0;j<nok;j++){
		    ///Recorremos todas las particulas que ya esten metidas en la lista
		    if(good){
		        r2=distancia2(bx, part, i, j);
                if(r2<1.1*0.5*(part[i].sigma+part[j].sigma)){
                    ///Hay solapamiento entre las particulas i y j. No aceptamos la particula i
                    good=false;
                    i--;
                }
                else
                {
                    if(r2<daux) daux=r2;
                }
		    }
		}
		if(good){

            if( dmin>daux) dmin=daux;
            ///Si pasa por aqui la particula no se solapa con ninguna otra --> se acepta
            nok++;
		}
	}
	cout<<"Particulas colocadas en la caja: (dmin: "<<dmin<<")"<<endl;

	r2=distancia2(bx, part, 160, 294);
	cout<<r2<<endl;

}


double anint(double x, double lx){
    double res=x;
    if (x>0.5*lx){
        res-=lx;
    }
    else
    {
        if (x<-0.5*lx){
            res+=lx;
        }
    }

    return res;
}


double distancia2(BOX bx, PART *part, int i, int j){
    double dis2=0;
    double drx,dry,drz;

    drx=part[i].x-part[j].x;
    while(drx>0.5*bx.lx)
        drx-=bx.lx;
    while(drx<-0.5*bx.lx)
        drx+=bx.lx;
    dry=part[i].y-part[j].y;
    while(dry>0.5*bx.ly)
        dry-=bx.ly;
    while(dry<-0.5*bx.ly)
        dry+=bx.ly;
    drz=part[i].z-part[j].z;
    #ifndef CONFINED
    while(drz>0.5*bx.lz)
        drz-=bx.lz;
    while(drz<-0.5*bx.lz)
        drz+=bx.lz;
    #endif

    double drx2=drx*drx;
    double dry2=dry*dry;
    double drz2=drz*drz;
    dis2=drx2+dry2+drz2;
    return  sqrt(dis2);
}
